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WIND- DRIVEN CURRENTS IN A SHALLOW LAKE OR SEA 


by Marvin Goldstein and Richard T Gedney 
Lewis Research Center 

ABSTRACT 

For shallow lakes and seas such as the great lakes (especially 
Lake Erie) where the depth is not much greater than the Ekman depth 
the usual Ekman dynamics cannot be used to predict the wind driven 
currents. The necessary extension to include shallow bodies of water, 
gi,yen by Welander, leads to a partial differential equation for the sur- 
face displacement which in turn determines all other flow quantities. 

A technique for obtaining exact analytical solutions to Welander's equa- 
tion for bodies of water with a large class of bottom topographies which 
may or may not contain islands is given. It involves applying conformal 
mapping methods to an extension of Welander's equation into the com- 
plex plane When the wind stress is constant (which is the usual assump- 
tion for lakes) the method leads to general solutions which hold for bodies 
of water of arbitrary shape (the shape appears in the solutions through 
a set of constants which are the coefficients in the Laurent expansion of 
a mapping of the particular lake geometry) The method is applied to 
an elliptically shaped lake and a circular lake containing an eccentrically 
located circular island 
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1, INTRODUCTION 

'! 

A main source of motion in the oceans andilakes is the wind. When 
the wind acts on the surface of a deep body of water, it sets up a circu- 
lation pattern which consists of a top and bottom friction (boundary) 
layer separated by a geostrophic (inviscid) core. In Ekman dynamics, 
the bottom stress is assumed proportional to the geostrophic velocity. 
However, Ekman dynamics is only valid when the thickness of the fric- 
tion (or Ekman) layer is negligible by comparison with the depth of the 
lake , The application of Ekman dynamics to any of the Great Lakes is 
questionable since they all have shallow regions of considerable extent. 
In Lake Erie, for moderate wind speeds, the thickness of the friction 
layer is comparable to the depth over much of the lake and Ekman 
dynamics can very definitely not be applied, 0 

The necessary extension of the Ekman analysis to the case of a 
homogeneous shallow lake has been given by Welander (1957) and that 
theory is used here for the case of steady winds 

t 

Gedney (1971) numerically calculated the wind-driven currents in 
Lake Erie using Welander *s formulation. Comparison of these calcula- 
tions with measurements was favorable and showed that the shallow 
sea formulation does indeed yield fairly accurate three-dimensional 
velocity profiles as a function of horizontal position and depth in the 
lake. The variation of bottom depth and boundary geometry must be 
included in any shallow lake analysis because the circulation depends 
very strongly on these factors. 
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Here a technique is given for obtaining exact analytical solutions to 
Welander’s equation for bodies of water which have bottom topographies 
of a general type and which may or may not contain an island. The shal- 
low lake model consists of a general second- order linear elliptic partial 
differential equation in two independent variables with coefficients which 
depend on the local bottom topography. This equation is transformed 
using conformal mapping methods into an equation whose coefficients 
are functions of only one of its two independent variables. By taking a 
finite Fourier transform the problem is reduced to the problem of solving 
a second- order ordinary differential equation. The solution is then repre- 
sented by an infinite series with each term a product of a function of one 
of the independent variables and a function of the other , When the wind 
stress is constant the method leads to general solutions which hold for 
bodies of water with shore lines of arbitrary shape. The shape appears 
in the solution through a set of constants which are the coefficients in 
the Laurent expansion of a mapping of the particular lake geometry 

2 , FORMULATION 

In the present analysis the basic approximations are that the water 
density is constant, the vertical eddy viscosity is independent of depth 
but dependent on wind velocity, the pressure is hydrostatic, and the 
lateral friction and nonlinear acceleration terms can be neglected. The 
above assumptions reduce the momentum equations to two equations con- 
taining the horizontal velocities and the surface slope as unknowns. The 
effect of the Coriolis force is included but its variation with latitude is 
not. taken into account. The appropriate boundary conditions for these 
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equations are a no- slip condition at the lake bottom and a specified shear 
stress (due to the wind) at the air- sea interface. These equations and 
boundary conditions can be solved analytically for the velocity 



where | is the vertical coordinate measured upward from the lake sur- 
face, £ is the surface displacement of the sea or lake above its zero- 
flow level, x and y are the horizontal, locally Cartesian coordinates, 
h is the depth of the body of water, d = 1^2 v/f c is the Ekman friction 
depth, v is the coefficient of vertical eddy diffusivity and f c is the 
Coriolis parameter which introduces the effect of the Earth's rotation. 

In addition U and V are the horizontal components of velocity in the 
x and y directions, respectively, g is the acceleration due to gravity 
and t w = + i where r™ and are, respectively, the x and 

y components of the wind stress divided by the density of water. 

By vertically integrating the horizontal velocity over the vertical 
direction and by using the continuity equation one can obtain a single 


equation for the surface displacement, which is 

so .so kin h\ i o __ h\l so ^ a ( h\ 


’xx b yy 




( 1 ) 
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where the subscripts x and y denote partial derivatives with respect 
to these variables and 



Equation (1) is a slightly rearranged form of that given by Welander 
(1957). For details of the derivation of the above equation and a discus- 
sion of the approximations involved in the derivation, see Gedney (1971). 
For instance Gedney has shown for Lake Erie that the approximations 
used, except for the approximation of a constant eddy viscosity, induce 
at the most a small (of the order of 10" *) error in the calculations. 

The x and y components of the total volume flow Q-^ and Qg 
respectively are related to the surface displacement by the relation 




gd 27T L 
where Q = Q 


ME + iF)« x + i? y ) 


( 2 ) 


1 + iQ 2 = J^ h 


(U + iV) d£ and 




The boundary conditions for equation (1) are obtained by specifying 
the volume flow normal to the coasts. Thus if n^ and ng are the x 
and y components, respectively, of the outward drawn unit normal n 
to the boundary, then for any closed body of water the boundary condi- 
tions for equation (1) are obtained by substituting equation (2) into the 
equation 

Qjn^ + Q2n2 — 0 (3) 

3. GENERAL SOLUTION 


We shall now suppose that there exists a nonconstant harmonic func- 
tion u of x and y and an arbitrary function H of u only such that 
the depth distribution h can be expressed in the form 


h(x, y)/d =— H 

2 7T 


3 u(x, y) 


(4) 


This is a fairly general functional form and it will be possible, for 
any one of a large number of lakes and seas, to choose the functions u 
and H in equation (4) so that the depth distribution is approximated 
fairly closely by a relation of this type. 

Now let v be the harmonic conjugate of u. Then the function 
w = u + iv (5) 

is an analytic function of the complex variable 
z = x + iy 


We are interested in obtaining solutions to equation (1) for a shallow sea 
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or lake which is either bounded by a single closed curve C such as that 
shown in figure 1(a) or a lake which contains a single island as shown in 
figure 2(a), Thus in the latter case we consider the doubly- corrected 
region where the outer boundary (shore line) is C and the inner (island) 
boundary is S. In either case we shall suppose that the bottom topogra- 
phy of the basin can be approximately described by a function of the 
forms (4) and that the depths along the shore line C in figure 1(a) and 
C and S in figure 2(a) are constants (which may in general be different). 
Equation (4) therefore shows that the harmonic function u must also be 
constant along C in figure 1(a) or along C and S in figure 2(a). With- 
out any loss of generality we can always redefine the functions H and u 
so that u - 0 on C. When an island is included in the lake we shall de- 
note the constant value of u on S by Ug. Again without loss of gener- 
ality, it is always possible to arrange matters so that Ug < 0. Now 
consider the analytic function w introduced in equation (5) and put 
T = e w 

Then T is an analytic function either within the simply- connected region 
of figure 1(a) or the multiply- connected region of figure 2(a).. Since 

it follows that J T | =1 for z on C and in the case where an island is 
i i U S 

present |T| =e = constant for z on S. Hence, when no island is 
present the mapping 
z — T 

transforms the simply connected region in figure 1(a) into the interior of 
the unit circle shown in figure 1(b) and in the case where a single island 
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is present the mapping 
z - T 

transforms the multiply connected region of figure 2(a) into the interior 

of the concentric annular region in the T- plane shown in figure 2(b), 

with the outer boundary C mapping into the unit circle and the inner 

Uo 

boundary S mapping into the inner circle which has radius R q = e 
As long as the boundaries of the regions shown in figures 1 and 2 
are sufficiently smooth, the Reimann mapping theorems for simply and 
doubly connected regions guarantee that the mappings z — T described 
above will always exist. 

We now introduce polar coordinates into the T- plane by the relation 


Since u is constant on the curve C and also when an island is present 
on the curve S, the boundary condition (3) becomes 

Q - 0 for zeC (6) 

\dw / 

and when the island is present 

YfL Ql— | = 0 for zeS (7) 

\dw/ 


Upon substituting the relation (4) into equations (1) and (2), then substi- 
tuting equation (2) into equations (6) and (7) and introducing r and 6 as 
the new independent variables we get 


r(rS r ) r + Z 0e + a(r)r? r 


+ b(r)C« =— 'Tk — 
y e(r) 0T 



( 8 ) 
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(for r = 1 ) 
f and for r = R Q i 
^island is present 


- IT — 9 ^ TT 


a + ib = i __ [sinh H(r ) - i sin H(r )] 2 r dH 

[cosh,H(r) + cos H(r)][sinh H(r) - sin H(r)] dr 

c 5 sinh H(r) - sin H(r) 1 

cosh H(r) + cos H(r) ( 

f _ sinh H(r) + sin H(r) 

cosh H(r) + cos H(r) v 

and the partial derivative with respect to T is taken with T* held con- 
stant. We have now shown that the surface displacement of the sea or 
lake can be found by solving equation ( 8 ) either in the unit circle of fig- 
ure 1 (b) if no island is present or in the concentric annular region of 

•> 

figure 2(b) if an island is present subject to the boundary condition (9). 

It now follows from equations ( 8 ) and (9) that if the complex function 
17 is the solution to the boundary value problem 

, . o„2 * 


:N r + + a<r)n, r + b(r) ^ = ^(Sy 


f(r) 


oT* / dz\* 


e(r) \dX 


(for r = 1 ) 
and for r = R q ' 
jf island is present 


- n ^ 9 ^ tt 


Then the surface displacement £ is determined by 

c=i6.d (is) 

The function r\ must certainly be periodic. Hence upon taking the 
finite Fourier transforms of equations (11) and (12) and substituting into 
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equation (13) we find that 

CO 

c ^o n (r)eine , (i4) 

n =- 00 

r(rft/)’ + ar£2^ - (n^ - inb)fl T n for n = 0, ±1, ±2, . . . (15) 

e 

(for r = 1) 


r «;- inia n =|y n { 

e e 




' and for r = R q \ 
if island is present/ \ 


where 


I =— (D T n — 
n 2niKj dT 


dT 


for n = 0, ±1, ±2,. . . 

(16) 

(17) 


o 



And the contour C Q is to be taken as any contour within the unit circle 
enclosing the origin for the case where no island is present and any con- 
tour within the annular region of figure 2(b) which encloses the origin 
for the case where an island is present. 

We have therefore reduced the problem of finding the surface dis- 
placement to that of solving the ordinary differential equation (15) sub- 
ject to the boundary condition (or conditions) (16). In the case where no 
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island is present the second boundary condition is replaced by a bound- 
edness condition on ft. at r = 0. 

4. CONSTANT WIND STRESS 

For lakes and to a lesser extent for shallow seas the case of prin- 
cipal interest for steady state winds is that for which the wind stress is 
constant. The reason for this is that the uncertainties in the wind veloc- 
ity measurements are of the order of the variation of this velocity over the 
body of water. In comparing measurements and calculations Gedney (1971) 
found this to be true for Lake Erie. In this case a is independent of 9 
and it follows from the Cauchy- Gour sat theorem that equations (18) and 
(19) become 

r =— r— (20) 

r n dr 

y n =— a (21) 

r 11 

It is easy to see from equations (20) and (21) that the boundary value 
problems (15) and (16) for ft n is independent of the actual shape of the 
body of water and depends only on the bottom topography. Hence the 
ft n can be determined once and for all for an arbitrary lake. The effect 
of the lake shape enters the problem in two ways. First equation (14) 
determines the surface displacement in terms of the parametric vari- 
ables r and 8. The mapping 
z - T 

which maps the particular shaped lake into the T- plane must then be 
used to relate r and 9 to the physical coordinate and therefore allows 
£ to be expressed as a function of x and y. 
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The effect of geometry also enters through sequences of constants 
{ I } . In fact it is easy to see that I is simply the coefficient of 
T~( n+ 1) in the Laurent series expansion of dz/dT about T = 0. This 
shows that for the case where an island is not present no positive power 
of n appear in the series (14). 

Thus, an exact solution for the sea- level elevation (and since all 
other physical quantities can be expressed in terms of this, an exact 
solution to the complete flow problem) can be obtained once the ordinary 
differential equation (14) with the boundary conditions (15) has been 
solved. 


5. SOLUTIONS WHEN THE DEPTH IS GREATER THAN 
ONE HALF THE EKMAN FRICTION THICKNESS 
It is shown by Goldstein, et al (1970) that for h/d s. 1/2 A, B, E, 
F, and a can be replaced with only a small loss of accuracy by 
A = 0 


B = 1 


E = 1 - 


27Th 

d 


F = 1 


r 2rtr w 2e -(l+i)H/2 
gd l 

Figure 3 shows how closely A, B, E, and F approximate these values 
for h/d s. 1/2. This approximation holds, for example, over the major 
portion of the great lakes. With these approximations equations (15) and 
(16) become 
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r(r£V )’ - (n - inrH’)12 n = -±- a' 

r 11 


( 23 ) 


rf2^ + in(H - 1) S2 n = — cr 

r 11 


(for r = 1) 

/and for r = R Q if\ 


(24) 


^island is present J 
In order to illustrate the method we shall work out specific examples 
using only the approximate equations (22) through (24). 

(a) Elliptic Lake Without An Island 
First consider the case of a lake with no island which has an approx- 


imately parabolic bottom topography given by 
H(r) = H q + 5(1 - r 2 ) 

where H Q + 6 is the constant depth at the shore. 


(25) 


Since we are considering the case where no island is present, I n is 


Hence only 12 n with negative values 


equal to zero for n = 0, 1, 2, • 
of n occur in the solution (14). When equation (25) is substituted into 
equation (22) we obtain an equation which is easily transformed into 
Bessel's equation. The solution to this equation which satisfies the 
boundary conditions (24) and which is bounded at r = 0 is given by 

>1 

-i-A 
^n 


7T 


12 =- A 

n 2 n 



— sC(r;p)a(p) dp for n = - 1, -2, ••• 


0 


n 


where 


<s^(r;p) s ( 


Y n<V> - "nVV*] J n+l<V> for » s r 
_ Y n+l<V> - "n J n+ l(V>] J n (X n r) for r s P 
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= 


n|l + i(H 0 


1) 


W' VWV 


' ■[■ 

A n =^/-2in5 


1 + i(H 0 - 1) 


J n^n^ " X n J n+l ( V 


and J and Y are Bessel functions, 
n n 

This solution can be applied to various shaped bodies of water simply 
by finding the conformal mapping z - T of the region occupied by the 
body of water in the physical plane (see fig. 1(c)) into the unit circle. 

Thus for an elliptically shaped lake of semi-major axis Lp and semi- 
minor axis L 2 the appropriate mapping which takes the center of the 
lake into the center of the unit circle (this causes the deepest portion of 
the lake to be located at its center) is 

T ='\/k sn| 


— sin*" 1 ( — 

n Vi/ 


, k 


where K is the complete elliptic integral of the first kind of modulus k, 
sn is the sine amplitude function 


I =-\ / juj ^ - ijg 

and k is determined from 


L 1 - L 2 \ 2 
L 1 + L 2/ 

The transformation T(z) could have just as easily been constructed so the 
deepest portion of the lake occurs at any point other than at the center of 
the lake. This is accomplished by performing an additional transforma- 
tion of the unit circle into itself with an arbitrary point mapping into the 
center and then forming the composite transformation. 


s -77K’/K _/ 



- 15 - 


(b) Circular Lake With An Eccentrically Located Circular Island 
Next consider the case of a lake with a single island which has a 
bottom topography given by 

H(r) - H +— — In r (26) 

0 In R 

o 


When this expression is substituted into equation (23) we obtain an equa- 
tion with constant coefficients. The solution to this equation subject to 
the boundary conditions (24) is 



for n = ±1, ±2, ±3, • • • 


where 

w n ( r ; 0 = [^ n - in(C - l)]r n + [A n + in(JJ - 1)] r n 
K n s K- in < H o ' iDOn + i"( H o + 5 ‘ ^ 


- I>„ + i"(H 0 - DJNi - ln(H o + 6 - R o " 

2n(2 + ia) -p+i)/2]ci! _ 1 


2n + (1 + i) (n + 


1 5 

2 In R 
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This solution can be applied to various shaped bodies of water and islands 
simply by finding the conformal mapping z — T of the doubly connected 
region occupied by the body of water in the physical plane into the concen- 
tric annular region in figure 2(b)- For example, consider a circular body 
of water with a radius R containing a circular island of radius R 

s 

whose center is located at some arbitrary point within the body of water. 
The configuration is shown in figure 2(c). In this case the appropriate 
mapping is 


T - z ~ aR 
az - R 


(28) 


where 



and the radius R Q of the inner circle in the T-plane is 



Upon substituting equation (28) into equation (17) and using the method of 
residues to evaluate the contour integral we find that 
CO for n = 0, 1, 2, • • • 




v R(a 2 - l)a n ” 1 


for n = - 1, -2, • • • 


( 29 ) 
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Substituting equations (27) and (29) into equation (14) determines the 
surface displacement. 

6. RESULTS 

A method for obtaining the analytical solution to equation (1) for the 
surface displacement in a lake or shallow sea has been given in terms 
of the power series 

? I* n n ( r) e in0 (14) 

n-_-oo 

The variables r and 0 are related to physical coordinates x and y 

of the lake through a conformal mapping function which depends only on 

the shape of the lake. The coefficient I* is determined solely by this 

mapping function. The bottom topography h(x, y) has been chosen in 

such a way that it depends only on r in the transformed r, 0-plane, 

In general the effect of the lake bottom topography as well as the lake 

boundary shape enters through the functions fi n (r) which are determined 

from the ordinary differentiation equation (15). However, for a constant 

wind stress O n (r) depends only on the bottom topography. 

Results will be given here for the case when the bottom depth is 

greater than one half the Ekman friction thickness and the wind stress is 

a constant. Under these restrictions the functions S2 (r) are determined 

n 

from the ordinary differential equation (23). Details of the determina- 
tion of the mapping function T(z), S2 n (r), and I* for each case are 
given in section 5. All solutions depend upon the wind shear stress, 
t w , the bottom depth, h(x, y), the friction thickness, d, and the lake 
boundary geometry. 
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The first example concerns an elliptically shaped body of water 
without an island and with the approximately paraboloidal bottom topogra- 
phy given by equation (25). The surface displacement which results for 
a depth variation which increases from a maximum of h/d = 9/2n at 
the center to a minimum of h/d = 5/2ir at the mainland shore is shown 
in figure 4. Here the wind direction is parallel to the x-axis and the 
ratio of semi minor to the semi major axis of the ellipse is 2/3. As is 
well known, the surface displacement for a constant depth basin with 
uniform wind stress is a plane inclined to the horizontal. The deviation 
of the surface displacement from a plane is then the effect of the varia- 
tion in bottom depth. The s- shape of the surface across the wind direc- 
tion in figure 4 will be found in most lakes because the dished- out type 
bottom topography used in this example occurs naturally. 

The second example considered is a circular lake containing an 
eccentrically located circular island with the logarithmic bottom topog- 
raphy given by equation (26). The results shown in figure 5 are for an 
island diameter of (x a - x^) = g R and an eccentricity of (x a + x^)/2R = 

0. 5. The bottom topography varies from a value of h/d = 1. 0 at the 
mainland shore to a value of h/d = 0. 5 at the island as shown in fig- 
ure 5(a). In an actual basin without an island, the depth is smallest at 
the mainland shore and largest at some central interior point. Islands 

i 

create local areas in the interior of the basin where the depth becomes 
shallow. This effect is modeled by the topography shown in figure 5(a). 

Figure 5(b) shows the surface displacement for the circular lake to 
be more of a plane than that shown in the previous case. The reason 
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for this is that the circular lake bottom topography is much flatter than 
the one used in the elliptic lake case. Although the deviations of the 
surface displacement from an inclined plane are small, they have major 
effects on the local velocities. 

The horizontal volume flow stream function is shown in figure 5(c). 
Here the stream function x p is related to the volume flows and Qg 
by 

Q x Q 2 =-^t 

1 dy z 3x 

The stream function plot consists of two gyres; the gyre on the right is 
rotating clockwise and the one to the left of the wind is rotating counter- 
clockwise. The two gyres are separated by a dividing streamline which 
in this case has a value of zero. The value of x}/ on both the island and 
mainland shores is zero and the zero streamline in the interior has been 
twisted relative to the wind in a clockwise direction. 

As shown by the xj; plots obtained by Gedney (1971), the dividing 
streamline for a dish shaped basin whose depth increases monotonically 
from the shore to some interior point always runs through the maximum 
depth point. If the depth decreases monotonically from the shore to an 
interior point, the dividing streamline would rim through the minimum 
depth point. In the case of the island shown in figure 5(c), the zero 
streamline runs through the minimum depth point interior to the mainland 
shore which happens to be the island boundary. If the island was not 
placed at the minimum depth point the value of xf; on the island boundary 
would be different from zero since the path of the streamline dividing the 
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two gyres is determined primarily by the extremes in the bottom topog- 
raphy. 

The horizontal velocities at |/d = 0, -0.125, -0.250, -0.375, and 
-0. 500 are shown in figures 5(d) through 5(h). The magnitude of each 
velocity vector plotted can be determined from the scale included on the 
plot. The origin of the velocity vector is the position at which velocity 
is actually occurring. The velocities at the surface are skewed to the 
right of the wind due mainly to the Coriolis force. Local perturbations 
in the surface velocities due to the volume flow gyres are evident. The 
acceleration of the flow around the island and the deceleration near 
x = -1 and +1 can be directly attributed to the volume flow gyres shown 
in figure 5(c). At |/d = -0. 125 mass is still being transported in the 
direction of the wind but an increasing amount is also being transported 
to the right of the wind. The flow acceleration around the island is still 
evident. At |/d = -0. 25 the flow pattern is very similar to that of the 
integrated volume flow shown in figure 5(c). Return flow opposite to the 
wind direction is now occurring. The flow near the island is predomi- 
nately tangent to it. At |/d = -0. 375 and -0. 50 there is much return 
flow opposite in direction to the top surface layer flow. The flow below 
|/d = -0. 50 is very similar in direction to but smaller in magnitude 
than that shown at |/d = -0. 50. This general flow pattern where mass in 
the top surface layer is transported along and to the right of the wind direc- 
tion and returned in the opposite direction in the bottom layer is the domi- 
nant pattern in shallow lakes. The variations in this general pattern are 
due to the particular bottom topography. For this particular case where 
the bottom slopes upwards at the island an acceleration of flow occurs 


around the island. 
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(d) Horizontal surface velocities. 
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(f) Horizontal velocities at f/d= -0.250. 


(g) Horizontal velocities at £/d« - 0 . 375 . 



( h> Horizontal velocities at £7d» -0.500. 
Figure 5. - Concluded. 
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